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Theoretical Models of the Galactic Bulge 


Juntai Shen & Zhao-Yu Li 


Abstract Near infrared images from the COBE satellite presented the first clear ev¬ 
idence that our Milky Way galaxy contains a boxy shaped bulge. Recent years have 
witnessed a gradual paradigm shift in the formation and evolution of the Galactic 
bulge. Bulges were commonly believed to form in the dynamical violence of galaxy 
mergers. However, it has become increasingly clear that the main body of the Milky 
Way bulge is not a classical bulge made by previous major mergers, instead it ap¬ 
pears to be a bar seen somewhat end-on. The Milky Way bar can form naturally from 
a precursor disk and thicken vertically by the internal firehose/buckling instability, 
giving rise to the boxy appearance. This picture is supported by many lines of evi¬ 
dence, including the asymmetric parallelogram shape, the strong cylindrical rotation 
(i.e., nearly constant rotation regardless of the height above the disk plane), the ex¬ 
istence of an intriguing X-shaped structure in the bulge, and perhaps the metallicity 
gradients. We review the major theoretical models and techniques to understand the 
Milky Way bulge. Despite the progresses in recent theoretical attempts, a complete 
bulge formation model that explains the full kinematics and metallicity distribution 
is still not fully understood. Upcoming large surveys are expected to shed new light 
on the formation history of the Galactic bulge. 


1 A brief overview on the properties of the Galactic bulge 


Most spiral galaxies consist of three main components, an invisible dark matter halo, 
an embedded flat disk, and a central bulge. The Milky Way is no exception. The 
Milky Way bulge comprises about 15% of the t otal luminosity, and its stellar mass 
is about 1.2 — 1.6 X IO'^’Mq ( Portail et al.L 2015h . Galactic bulges contain crucial in¬ 
formation about the galaxy formation and evolution. Major mergers between galax- 


Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Obser- 
vatoiy, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China, e-mail: 
jshen,lizy@shao.ac.cn 


1 







2 


Shen & Li 


ies generally create a significant classical bulge that is similar to ellipticals in many 
aspects (see Chapter 6.1 of this book), whereas the long t erm internal secular evo¬ 
lution in the disk galaxy tends to build up a pseudobulge ( Kormendv & Kennicutd 
20041) . Understanding the structure of our Milky Way bulge is nontrivial, mostly 
because of our location in the disk plane and the severe dust extinction in the op¬ 
tical band. One the other hand, a huge advantage being inside the Milky Way is 
the power to resolve and observe individual stars. Here we briefly summarize the 
structure, chemical composition, age, and kinematics of the Galactic bulge. More 
d etailed reviews o n the observ a tional properti es of the Galactic bulge can be found 


in 


Zoccali ( 2010l) : Rich ( 2013 ): Origlia ( 2014 ) and Chapter 4.1 of this book. 


Structure. The pre sence of a triaxia l structure in the inner Galaxy was first hinted 
from gas kinematics ( de Vaucouleurs . 1964 : Binnev et all 199ll: Burton & Liszt . 
19931). The near infrared images from the COBE satellite revealed clearly that 
the Mil ky Way conta i ns an asymmetric parallelogram-shaped boxy bulge in the 


center (IWeiland et al.L 119941) . The asymmetry may be explained by a tilted bar; 
the near end of the bar is closer to us than the far side, consequently it ap- 


txx>_^ L^LXX XC X..XWCWX CXC LXXLAXX LXX>-^ X LAX XL LX^ 

pears to be bigger than the other side ((Blitz & Spergeli Il99lh . Measurements of 
the three-dimensional density distribution of the Galactic bulge using various stel¬ 
lar tracers give a triaxial bar with the semi-major axis ^ 3 — 4 kpc and tilted by 
^ 20° — 30° between the Sun-Galactic Cent er (GC) line (i.e., StaneketakL 1997 
Bissantz & Gerhard. 20021: Rattenburv et ak. 2007bt Robin et al. . 2012 : Cao et al 


2013t Wegg & Gerhardl 2013t Pietrukowicz et al. . 2014). This main triaxial bar is 
sometimes dubbed as “the (boxy) bulge bar” (i.e., boxy pseudobulge in the notation 
of this book), whose properties still differ among different research groups. 

The existence of other bar structures than the main bulge bar is still under active 
debate. Based on stars coun ts from the GLIMPSE ( Galactic Legacy Mid-Plane Sur¬ 
vey Extraordinaire) survey, Beniamin et al. ( 20051) argued for another planar long 
bar passing throu gh the GC with half-length 4.4 kpc tilted by ^ 45° to the Sun- 
GC line (see also Cabrera-Lavers et al. 2007h . If this long bar is confirmed, then 
its co-existence with the similarly-sized bulge bar is dynamically puzzling, as their 
mutual torque tends to ali gn the two bars on a short time scale. Also both observa¬ 
tions of external galaxies ( Erwin & Sparkel 20021 2003 ) and simulations of long- 
lived double barred galaxies ( Debattista & Shen . 20071: ^en & Debattistal 2009h 
have shown that the size ratio of two bars is genera l ly abo ut 0.1 — 0.2. Using nu¬ 
merical simulations, Martinez-Valpuesta & Gerhardl ( 2011 ) showed that the obser¬ 
vational signatures of the planar long bar may actually be reproduced with a single 
bar structure. They suggested that the long bar may correspond to the leading ends 
of the bulge bar in interaction with the adjacent spiral arm heads. Nishivama et al.l 
(2005); Gonzalez et al.l ( 201 lat) also suggested the existence of a possible secondary 
inner bar, as the slope of the longitudinal magnitude peak profile of red clumps flat¬ 
tens at |/| ^ 4°, deviating from the main bar . With the same single barred A-body 
model, Gerhard & Martinez-Valpuestal ( 2012l) showed that such a slope change may 
be caused by a transition from highly elongated to nearly axisymmetric isodensity 
contours in the inner boxy bulge. 
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Chemical composition and stellar age. The chemical composition and the stel¬ 
lar age are crucial parameters to constrain galaxy formation history and to estab¬ 
lish the connection to stellar populations in other structural c omponents. The bulk 
of bu l ge stars is old with a wide range of metal abundances (IMcWilliam & Richl 


19941 Zoccali^^ n 20081). including some of the o ldest stars in the Milky Way 
(e.g. Howes et ah . 20141 Schlaufman & Casev . 2014). The metallicity distribution 
in the bulge displa ys a vertical gradient both on and off the minor axis away from 


the Galactic plane ( Zoccali_et^ , 200^ Gonzalez et al. , 2011b ; .lohnson et al. , 2011 


2 OI 2 I 2013h . while Rich. Origlia & Valenti ( 2012 ) found no major ve rtical abun¬ 


dance gradient close to the disk plane {b < 4°). iNess et al.l (l2013al) found that 


stars with [Fe/H] > —0.5 are part of the boxy bar/bulge, and the metal-poor stars 
are likely associated with thick disk. As shown in Fig. [1] although the rotation 
curves are similar for three different metallicity bins, the metal-poor population 
(—0.5 < [ Fe/H] < — 1 .0) clea rly has higher velocity dispersion and spheroidal kine¬ 
matics. In iNess et al. ( 2013a ). stars with [Fe/H] ^ 4-0.15 are more prominent close 
to the plane than the metal-poor stars, appearing as a vertical abundance gradient 
of the bulge. Most bulge metal-poor stars show enhanced alpha elements com¬ 
pared to both thin and thick di sk stars ( .Johnson et al.L 2011; Gonzalez et Iil l2011bl: 
Rich. Origlia & ValentiL 2012 ). suggesting a rapid formation timescale with respect 


of both disk components. 

Compared to metallicity, age determination is more difficult and imprecise (see 
the review by ISoderblom 2010 and referen c es therein). The bulk of bulge stars is 


old ( ^ 10 Gyr) (e.g., lOrtolani et all Il99-5t iLecureur et all l2007t IClarkson et al 


20081). However, intermediate-age metal-rich stars are also detected in the bulge 


region, with t he exact relative fraction still under debate (IClarkson et all 1201 1 


Bensbv et al. . 2011 , 2012 : Ness et all 2014 ). In addition, there is a nuclear disk 


of much younger stellar population with ongoing star formation in the central 
200 pc (sometimes termed “nuclea r bulge”), whose mass is about 1.5 x IO^Mq 
( L aunhardt. Zvlka & Mezg^ 2002 ). 

Kinematics. To study systematically the stellar kinematics, iRich et al.l (l2007h 
initiated the Bulge Radial Velocity Assay (BRAVA) project with M giants as tracers 
covering the whole Galactic bulge (~9000 stars). The BRAVA survey revealed clear 
cylindrical rotation of t he bulge, i.e., nearly constant rotation regardle s s of the height 


above the disk plane (iHoward et al.l l2008l l2009t iRich et al.L 120081 : iKunder et al 


2012h . BRAVA kinematics also put th e Galact i c bulg e close to the oblate isotropic 
rotator line in Vmax/o’— e diagram ( Binnevl 1978 ). distinct clearly from a hot 
slowly-rotating system like the Milky Way halo supported by velocity dispersion. 
The BRAVA radial velocity distribution was well reproduced by a self-consistent 
V-body model of a pure-disk Galaxy by Shen et al. (2010), with no need for a sig- 
nihcant classical bulge (see Section|2|i. 

The Abundance and Radial velocity Galactic Origins Survey (ARGOS) obtained 
radial velocities and stellar parameters for 28000 stars in the bulge and inner disk 
of the Milky Way galaxy across latitudes of b = —5° and —10° dFreeman et al.L 
2013h. The cylindrical rotation of the bulge was also conhrmed in the ARGOS 


data (iNess et al.l l2013al) . They found a kinematically distinct metal-poor popula- 
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Fig. 1 Rotation and velocity dispers ion profiles in the ARGOS observations towards the Galactic 
bulge fields from lNess et alJ ll2013ah . The three columns correspond to three different metallicity 
bins, decreasing from le ft to right. Different symbols represent stars in different fields. Reproduced 
from \Ness et g/J il207jd) . 


tion ([Fe/H] < —1.0), which may be related to the thick disk or halo. In the com¬ 
missioning observation of APO Galactic Evolution Experiment (APOGEE) towards 
the Galactic bulge region, a high Galactocentric velocity (Vgsr -1-200 km s^*) 
and cold ( g ^ 30 km sy*) stre am was re ported and sug gested as the bar support¬ 
ing orbits ( Nidever et al. , 2012h . However. Iki et al.l ( 2014 1 suggested that typical bar 
models do not generate a distinct cold high velocity stream observed from the so¬ 
lar perspective. A smooth high velocity shoulder instead does exist in many bulge 
fields; it roughly corresponds to the tangential point between the line-of-sight and 
the bar-supported orbits as shown in the distance-velocity diagram. The cold high 
velocity stream, if confirmed, may not be necessarily associated with the bar, but 
could be due to other substructures. Recent Giraffe Inner Bulge Survey (GIBS) ob¬ 
served 24 Galactic bulge fields and confirmed the cylindrical rotation and the lack of 


a sign ificant cold high velocity peak in the radial velocity distribution (IZoccali et al 


2014h. Besides the radial veloc ity, the proper motion is also an important parameter 


(e.g., Rattenburv et al. . 2007 ah . Given the stellar distance, the transverse velocities 
can be estimated. In the Baade’s window, the velocity ellipsoid of metal-rich stars 
shows a vertex deviati on in the radial versus transv erse velocity, consistent with the 
bar supporting orbits ( Soto. Rich & Kuiikenl 2007h . 

Ba sed on the radial veloc ity and the [Ee/H] measurement of three minor axis 
fields, Babusiaux et al. ( 2010h identified two distinct population; the metal-rich pop¬ 
ulation with bar-like kinematics and the metal-poor pop ulation corresponding to 
an old spheroid or a thick disk (also see Hill et al. . 201 ih . The metal-rich popula¬ 
tion demonstrates smaller velocity dispersion and lower alpha-element enhancement 
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comp ared to the metal-poor population ( Johnsonet all |201 it lUttenthaler et al 


20121) . Across the bulge fields in ARGOS survey, the metal-poor population ([Fe/H] < 


— 1.0) was also found to be kinemat i cally di stinct with large velocity dispersion and 


non-cylindrical rotation (INess et al.Ll2013bl) . 


X-shaped structure. A recent development in the bulge structural study is the 
discovery of an intriguing X-shaped structure (Section (3). Recently, two groups 
independently reported the bimodal brightness distrib ution of the red clump (RC ) 
stars, which can be cons idered a good standard candle dStanek & Garnavi^ ^ 1998 ), 


in the Galactic bulge ( MeWilliam & Zoccalil 2010l hereafter MZIO; Nataf et al 


2010l) . MZIO suggested that the bimodality is hard to explain with a naive tilted 


bar since the line of sight crossing the bar can only r esult in stars with one dis¬ 
tance. One possibility speculated by iNataf et aP (l2010l) is that one RC population 
belongs to the bar and the other to the spheroidal component of the bulge. Another 
puzzling fact is that distances of the bright and faint RCs are roughly constant at 
different latitudes, which was hard to understand with a naive straight bar. They 
proposed that these observed evidences can be well explained with a vertical X- 
shaped structure in the bulge regio n. The existence of this particular structure was 
later verified bv ISaito et alJ (1201 ll) . They found that the X-shaped structure exists 
within (at least) |/| <2°, and has front-back symmetry. Around b = ±5°, two RCs 
start to merge, due to s evere dust extinction and foreground contaminati on (MZIO; 
Wegg & Gerhardl2013 ). From numerical simulations, as demonstrated in Li & Shen 


( 2012h an'd Ness et al.l ( 2012 l. the buckled bar naturally reproduces the observed X- 
shape properties in many aspects. The X-shape extends to about half the bar length 
with similar tilting angle as the bar. The observed north-south symmet ry of the X- 
shape indicates that it must haye formed at least a few billion years ago ( Li & ShenL 


20121 ). 


2 A fully evolutionary bar model as the basis of understanding 
the Galactic bulge 


Theoretical modeling of the Milky Way bulge made intense use of A-body simu¬ 
lations. The basis of a successful Galactic bulge model is a fully eyolutionary bar 
model that deyeloped naturally from the bar instability of a cold massiye disk. Here 
we describe a successful high-resolution bar model deyeloped in Shen et al. ( 2010[ 
hereafter S10), which was initially motiyated to match the BRAVA stellar kinematic 
data. A-body bar m odels to exp l ain the Galactic bulge were already attempted in 
early studies such as FuxI ( 1997 ): Seyenster et al. ( 19991) . but little stellar kinematic 
data were ayailable to constrain their models. 

SIO simulated the self-consistent formation of a bar that buckles naturally into 
a thickened state, then scaled that model to fit the BRAVA kinematic data on bulge 
rotation and random yelocities. BRAVA measured the stellar radial yelocities of 
M-type giant stars whose population membership in the bulge is well established. 
These giants proyide most of the 2 fim radiation whose box-shaped light distribution 
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motivates bar models. SIO used nearly 5,000 stellar radial velocities in two strips at 
latitude b = —A° and b = —8° and at longitude —10° < I < +10°, and a strip along 
the minor axi s (/ = 0°). The strong cylindrical rotation was found in the preliminary 
BRAVA dat a ( Howard_et_^ 2009ll consistent with an edge-on, bar-like pseudobulge 
( Kormendv 19931: Kormendv & Kennicutt . 20041) . although a precise fit of a bar 
model to the data was not available. The success prompted SIO to construct a fully 
evolutionary V-body model that can fit the radial velocity data of BRAVA. 


t= 900 Myr 



Fig. 2 Sh ortly after its formation, the bar becomes vulnerable to the vertical firehose/buckling 
instability l lToomrel . ll966l : lRaha et alill99lh . The figure illustrates the vigorous buckling instability 
that makes the initially thin bar bend out of the disk plane, reaching a considerable maximum 
distortion. After the buckling instability saturates on a short dynamical timescale (in a few hundred 
million years), the bar is greatly thickened in the vertical direction, giving rise to the boxy shape, 
(see Fig.|4). 


SIO used a cylindrical particle-mesh code ( Shen & Sellwood , 20041) to build fully 
self-consistent V-body galaxies. The code is highly optimized to study the evolution 
of disk galaxies, and it allows SIO to model the disk with at least 1 million particles 
to provide high particle resolution near the center where the density is high. They 
tried to construct the simplest self-consistent V-body models that fit the BRAVA 
data, avoiding contrived models with too many free parameters. Initially, they con¬ 
tained only an unbarred disk and a dark halo. The profile of the Galactic halo is 
poorly constrained observationally; SIO adopted a rigid pseudo-isothermal halo po¬ 
tential which gives a nearly flat initial rotation curve between 5 kpc and 20 kpc. A 
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Fig. 3 The horizontal line marks the pattern speed Qp of the quasi-steady bar in internal simulation 
unit with = G = = 1. Here Qp Ki 39 km/s/kpc in physical units. The solid line shows the 

curve of the circular a ngula r frequ ency Q, and the dashed lines mark Q ± k/2 at around t = 4.8 
Gyr. Reproduced from \Sheri il20i^) with permission. 


simple halo form allows SIO to run many simulations quickly, which greatly facili¬ 
tates a parameter search. A rigid halo also omits dynamical friction on the bar, but 
the central density of this cored halo is low enough so that friction will be very mild. 
Since the bulge is embedded well interior to the core radius of the halo, the exact 
profile of the dark halo at large radii is not critical. 

The process of bar formation and thickening in SIO is physically understood. 
First a bar develops self-consistently via the bar instability from the initially un¬ 
barred, thin disk. Bar formation enhances the radial streaming motions of disk par¬ 
ticles, so the radial velocity dispersion quickly grows much bigger than the verti¬ 
cal one. Consequently the disk buckles vertically out of the plane like a firehose 
.l^fl this i s the well known firehos e or buckling instability (e.g.. lToomre[ll965 


Combes et al.L 1990l : Raha et all 199 ih . It raises the vertical velocity dispersion and 


increases the bar’s thickness. This happens on a short dynamical timescale and sat¬ 
urates in a few hundred million years. The central part of the buckled bar is elevated 
well above the disk mid-plane and resembles the peanut morphology of many bulges 
including the one in our Galaxy. 

SIO found the one that best matches our BRAVA kinematic data after suitable 
mass scaling, out of a large set of V-body models. The barred disk evolved from a 
thin exponential disk that contains = 4.25 x about 55% of the total mass 

at the truncation radius (5 scale-lengths). The scale-length and scale-height of the 
initial disk are ~ 1.9 kpc and 0.2 kpc, respectively. The disk is rotationally supported 
and has a Toomre-Q of 1.2. The amplitude of the final bar is intermediate between 
the weakest and strongest bars observed in galaxies. The bar’s minor-to-major axial 
ratio is about 0.5 to 0.6, and its half-length is ~ 4 kpc. Fig. [3 shows the pattern 
speed of the bar {Qp w 39 km/s/kpc) and locations of the Lindblad resonances. The 


* It is puzzling why no galaxy has been caught in the process of violent buckling as shown in Fig.|3 
The saturation timescale of the buckling instability is very rapid (about a few hundred Myr), but it 
is not short enough to miss out the violent buckling phase if one can observe thousands of edge-on 
barred galaxies. Perhaps this implies that the buckling instability must have happened in the very 
early assembly stage of disk galaxies, which is much harder to observe. 
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bar properties in the best-fitting model are comparable to those obtained in other 
independent studies. The left panel of Fig. |4] shows face-on and side-on views of 
the projected density of the best-fitting model in SIO. A distinctly peanut shaped 
bulge is apparent in the edge-on projection. Fig. 0] (right panel) shows the surface 
brightness distribution in Galactic coordinates in solar perspective. Nearby disk stars 
dilute the peanut shape, but the bar still looks boxy. The asymmetry in the longitu¬ 
dinal direction can be understood as a perspective effect; the near end of the bar 
(at positive Galactic longitude) is closer to the Sun, so it appears bigger and taller 
than the far side. Both the boxy shape and the asymmetry are in good agreement 
with the morphology revealed by the near-infrared image from the COBE satellite 
dWeiland et al.Lll994 . 

Fig. |5]compares the best-fitting model kinematics in SIO (solid lines) with the 
mean velocity and velocity dispersion data from the BRAVA and other surveys. All 
velocities have been converted to Galactocentric values (the line-of-sight velocity 
that would be observed by a stationary observer at the Sun’s position). For the first 
time, this model is able to simultaneously match the mean velocities and velocity 
dispersions along two Galactic latitudes (—4° and —8°) and along the minor axis. 
The model comparis on with the complete BRAVA data release was also impressive 
( Kunder et all 2012 ). 

SIO also provided some constraints on the bar angle between the bar and the 
Sun-GC line. Using both the fit to the velocity profiles and the photometric asym¬ 
metry, they found that the overall best-fitting model has a bar a ngle of ^ 20° , which 
agrees reasonably well with other independent studies (e.g., StaneketakL 1997 : 


Rattenburv et al.. boOTbl: Cao et al., 20131: Wegg & Gerhardl 2013h. 

The best-fitting model in SIO contains no classical bulge component. SIO also 
tested whether or not a significant classical bulge is present, since it could have 
been spun up by the later formation of a bar, flattened thereby and made hard to 
detect. They found that including a classical bulge with >~ 15% of the disk mass 
considerably worsens the fit of the model to the data, even if the disk properties are 
accordingly re-adjusted. If the pre-existing classical bulge is overly massive, then it 
becomes increasingly ha rd to match both the mean velocity and veloc ity dispersions 
simultaneously (see also Saha, Martinez-Valnuesta & Gerhard 2012 ). 

The BRAVA kinematic observations show no sign that the Galaxy contains a sig¬ 
nificant merger-made, classical bulge. SIO demonstrated that the boxy pseudobulge 
is not a separate component of the Galaxy but rather is an edge-on bar. This re¬ 
sult also has important implications for galaxy formation. From a galaxy formation 
point of view, we live in a pure-disk galaxy. Our Galaxy is not unusual. In fact, gi¬ 
ant, pure-disk galaxies a re common in environments like our own that ar e far from 
rich clusters of galaxies ( Kormendv et al. . 201f)t Laurikainen et al. . 2014 ). Classical 
bulgeless, pure-disk galaxies still present an acute challenge to the current picture 
of galaxy formation in a universe dominated by cold dark matter; growing a giant 
galaxy via hierarchical clustering involves so many me rgers that it seems almos t 
impossible to avoid forming a substantial classical bulge dPeebles & Nussei . 20101) . 
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Fig. 4 Left three panels: face-on and side-on views of the surface density of our hest-fitting model 
as seen from far away. The Sun’s position 8.5 kpc from the Galactic center is marked along the 
+x axis. The Galaxy rotates clockwise as seen in the face-on projection. Right panels: The COBE 
DIRBE composite image of the Milky Way bulge and model’s surface brightness map in solar 
perspective. Both the model and observed image show the box-shaped, edge- on bar that appear 
bigger on its near side (positive longitude side). The left panel is reproduced from \Shen et al\ i20ld) 
with permission. 



Fig. 5 (top): Mean velocity and velocity dispersion profiles of the best-fitting model (black lines) 
compared to all available kinematic observations. The left two panels are for the Galactic latitude 
b = —4° strip; the middle two panels are for the b = —8° strip; and the right two panels are for 
the 1 = 0° minor axis. The black dia monds and their error bars are the B RAVA data; the green 
diamonds are for M-type giant stars ^Rangwda,^illiams_&^tmied|200^, and the red triangles 
are the data on red clump giant stars iRangwala. Williams & StanekLl200^ . This is the first time 
that a single dynamical model has been co mpared with data of such quality. The agreement is 
striking. Reproduced from \Shen et ai\ ilOldj) with permission. 
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2.1 The bulge structure may be younger than its stars 


A common misconception about the secularly evolved bar model is that it seems in¬ 
consistent with the old age of bulge stars that formed on a short timescale (~ 1 Gyr), 
as demonstrated by their a-element enhancement. It is important to make a distinc¬ 
tion between the assembly time of the bulge/bar and the age of the bulge stars; the 
Galactic bulge structure may be younger than its stars. The stars in our Galactic bar 
are older than most disk stars, but those star s could have formed over a sh ort period 
of time lo ng before the bar structure formed ( , 199^ Freem^ 2008h. Their old 


age (e .g.. lZoccali et al.Ll2003tlFulbright. McWilliam & Richl I2007t iClarkson et al 
20081) is therefore not an argument against the internal secular evolution model. 


2.2 Could the vertical metallicity gradient be produced in a simple 
bar model? 


Although the general bulge morphology and kinematics are well explained by 
this simple and self-consistent model, a potential difficulty of the simple model is 


how to explain the obse rved vertical metallicity gradient (e.g.. lZoccali et al.L 12008 


Gonzalez et all 201 Ibh . Intuitively one may expect that any pre-existing vertical 
metallicity gradient should be erased due to mixing in the violent buckling process. 
Thus the vertical metallicity gradient has been suggested as the strongest evidence 
for the existence of a classical bulge in the Milky Way. 

SIO proposed a plausible solution that an abundance gradient can still be pro¬ 
duced within the context of secular bar/bulge formation if some of the verti- 
cal thickening is produced by resonant heating of stars that scatter off the bar 
dPfenniger & NormanLll990r) . If the most metal-poor stars are also the oldest stars, 
then they have been scattered for the longest time and now reach the greatest heights 
away from t he disk plane. _ 

Recently Martinez-Valnuesta & GerhardI ( 2013 ) challenged the widespread be¬ 
lief that secularly evolved bar/bulge models cannot have metallicity gradients sim¬ 
ilar to those observed in our Galaxy. Using a similar boxy bulge/bar simulation 
as in SIO, they were able to successfully reproduce the observed vertical metallic- 
ity gradient and long itude-latitude metallicity map similar to that constructed by 
Gonzalez et alJ ( ^13l) . provided that the initial unbarred disk had a relatively steep 
radial metallicity gradient. One important assumption is that the pre-existing radial 
metallicity gradient is set up during the buildup of the precursor disk prior to bar 
formation. They proposed that if the Galactic bar/bulge formed rapidly from the 
precursor disk at early times, the violent relaxation may be incomplete during the 
bar and buckling instabilities, thus transforming radial pre-existing metallicity gra¬ 
dients to vertical gradients in the hnal boxy bulge. They further proposed that the 
range of bulge star metallicities at various latitudes may be used to constrain the 
radial gradient in the precursor disk. In their study, they tagged particles with some 
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metallicity and no chemical evolution was considered. If their result is confirmed 
in more detailed studies, then the vertical metallicity gradient is no longer a strong 
argument against the secularly-evolved bar/bulge model. 


2.3 Merits of the simple self-consistent bar/bulge model 


The main advantage of a simple self-consistent bar/bulge model such as the one 
in SIO lies in its simplicity; SIO tried to construct the simplest self-consistent N- 
body models that fit the BRAVA data, avoiding contrived models with too many 
free parameters. In addition, the SIO model is not just a static model but rather one 
that evolved naturally to this state from simple initial conditions. The bar is self- 
consistently developed from a massive cold precursor disk embedded in a cored 
halo. So it has relatively few parameter to tweak, unlike the more complicated 
chemo-dynamical models. 

Secondly, physical processes dictating the formation and evolution of the bar/bulge 
can be well understood. The formation of evolution of the bar is a natural outcome 
of two well-studied dynamical instabilities, namely the bar instability creating the 
bar and the subseque nt firehose/buckling instability that thickens the bar vertically 
(see SellwoodI 20141 for a comprehensive review on these instabilities). The buck¬ 
ling instability also helps the Milky Way to develop an intriguing X-shaped structure 
(Section|3ll. The different components of the early inner Galaxy (thick disk, old and 
younger thin disks) probably become trapped dynamically within the bulge structure 
kess et al.Lr2013bh . 

Thirdly and probably most importantly, despite its simplicity, the best-fitting 
model of SIO also ties together several isolated results (e. g., photometric bar angle, 
kinematic bar angle, successful kinematic fits in the context of observed cylindrical 
rotation, reasonable bar length and bar pattern speed, resulting constraint on a clas¬ 
sical bulge, the vertical metallicity gradient, and explaining the X-shaped structure) 
into a coherent picture, cemented by an A-body model in which the bar evolves nat¬ 
urally and without complicated fine-tuning. At current stage, it is still unrealistic to 
expect a full chemo-dynamical model to match all observed properties of the Galac¬ 
tic bulge. The simple bar model was designed to serve as a physically-motivated 
starting point, then one can gradually incorporate more complexities of the Milky 
Way bulge on top of it. 


3 The X-shaped structure in the Galactic bulge 


RCs are a good standard candle commonly used in the studies of the Galactic bulge. 
The apparent magnitude of RCs in high-latitude bulge fields shows a clearly bi- 
modal distance distribution, often termed “split r ed clump”. This b imodal distribu¬ 
tion was discovered in the 2MASS data (MZIO, Saito et al. 2011) and the OGLE 
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Fig. 6 a). Left panels: Distance distributions of particles in the SIO model in the Galactic bulge 
fields at a given longitude (top row) and latitude (bottom row); b). Right panels: Demonstration of 
the X-shape structure with the upper panel showing the side-on view of the bar in the SIO model 
and the lower panel showing the residual after fitting and subtracting the underlying smooth light 
contribution. The vertical X-shap ed structure is hig hlighted in this residual image. The length unit 
is Ri = 1.9 kpc. Reproduced from \Li & Shem\201^) with permission. 


data (Nataf et al. . 2010h . then confirmed also in the ARGOS sample ( Ness et al 
2012h . The RCs seem to be distributed in a vertically-extended X-shaped structure. 
This X-shape was initially puzzling; it seemed difficult to explain them using a tilted 
naive ellipsoidal bar. 

This X-shaped structure does not have a straight-forward explanation in classi¬ 
cal bulge formation scenarios, but it is a natural consequence of the bar buckling 


processes (see SectiorlUl if it is properly mode l ed ( Combes_et^ , 199^ Ruhaet^ 


1^91 ; Bureau&Ffeem^ , 1999t Athanassoulal 2005 : Martinez-Valpuesta. Shlosman & Heller 


2006). Li & ShetT(2012h analysed the same best-fitting Milky Way bar/bulge model 

in SIO, and demonstrated that it can qualitatively reproduce many observations 
of the X-shaped structure, such as double peaks in distance histograms (MZIO, 
Nataf et ^l2010h an d number density maps dsaito et ^ 2011 ). 


Li & ShenI (12012h found that an X-shaped structure is clearly discernible in the 


inner region of the side-on view of the SIO best-fitting bar/bulge model (top panel 
of Fig. IQr). They fitted and subtracted the underlying smooth component from the 
side-on bar model, then the X -shaped structure is highlighted in the bottom panel 
of Fig. [hjt dLi & ShenL 2012h . Fig. |6^ shows the distance distributions of particles 
towards the Galactic bulge, where double peak features similar to observations can 
be clearly identified. The bottom row shows that at a given latitude b the peak posi¬ 
tions are roughly constant. As the longitude decreases, the peak at larger distances 
becomes stronger. The top row shows that at a given longitude I the separation of 
the two peaks increases as the line of sight is further away from the Galactic plane. 
These results nicely agree with the observations of the X-shaped structure (MZIO). 
Along the bar major axis, the end-to-end separation between the inner two edges of 
the X-shaped structure is ^2 kpc. For the outer two edges, the end-to-end separation 
is ^4 kpc. The size of the X-shaped structure along the bar is estimated by averaging 
the two separations, which yields ^^3 kpc. It is worth noting that the value is much 
less than the bar’s full length (^8 kpc). Similarly, the end-to-end vertical separation 
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between the inner two edges of the X-shaped structure is '^1.2 kpc. For the outer 
two edges in the vertical direction, this separation is ^2.4 kpc. Therefore the verti¬ 
cal size of the X-shaped structure in the S10 model is '^l.S kpc. By s umming up the 
pixels with positive values in the X-shaped region, Li & Shenl ( 2012 ) estimated that 
the light fraction of this X-shaped structure relative to the whole boxy bulge region 
is about 7%. It is still unc ertain how much m ass is contained in the X-shape. More 
sophisticated analysis by IPortail et alJ (120 15h . based on the reconstructed volume 
density of the Galactic b ulge from Vista Variables in the Via Lactea survey (VVV) 
( Wegg & Gerhardl 2013h . suggests an off-cen tered X-shape enclosing about 20% of 
the bulge mass. Recently Nataf et al. ( 2015h studied the X-shape properties based 
on OGLE-III observa tions, and also found good agreement with the models in SIO 
and Ness_et^ (2012). 

Li & Shen ( 2012h also demonstrated that the X-shaped structure becomes nearly 


symmetric with respect to the disk plane about 2 Gyr after the buckling instability 
gradually saturates. The observed symmetry (MZIO) probably implies that the X- 
shaped structure in the Galactic bulge has been in existence for at least a few billion 
years. 


Based on the ARGOS sample. iNess et al.l (120121) found that the X-shaped struc¬ 
ture isjnainlyxomposedof the metal-rich stars r ather than the metal-poor ones (also 


see 


Uttenthaler et al. 201J^ Gihnoree^^ 201^. However, this conclusion is ques¬ 
tioned by Nataf. Cassisi & Athanassoulal ( 2014h who demonstrated that there may 


be a bias of metallicity on the RCs distance determination. In the future larger and 
unbiased samples are required to answer this question unambiguously. 

The existence of the X-shaped structure in our Milky Way provides extra evi¬ 
dence that the Galactic bulge is shaped mainly by internal disk dynamical instabil¬ 
ities instead of mergers, b ecause no other known physical processes can naturally 
develop such a structure. De Propris et alJ ( 2011 ) studied the radial velocity and 
abundances of bright and faint RCs at (/,h) = (0°,—8°), and found no significant 
dynamical or chemical differences. This may suggest that the two RCs indeed be¬ 
long to the same coherent dynamical structure, which can be naturally made in the 
formation of the bar/boxy bulge. 

Orbital structure studies are essential for understanding the properties of the X- 
shaped structure, which is the outcome of the collective buckling instability. The 
backbone orbits of a three-dimensional buckled bar are the xi tree, i.e., the x\ fam¬ 
ily plu s a tree of three-dimensional families bifurcating from it (IPfenniger & Friedli , 
199 ih . It is widely believed that the X-shaped structure may be supported by orbits 
trapped around the three-dimensional xi family (a lso known as banana orbits due 
to thek banana shape when viewed side-on) (e.g., Patsis., Skokos & Athanassoula , 
2002t Skokos. Patsis & Athanassoul4 20o3^ Patsis & Katsanika4 2014al). But this 


picture was recently challenged bv IPortail Wegg & Gerhardl (l2015l) w ho argued 
that th e X-shaped bulge is mainly supported by brezel-like orbits instead. Pin et al.l 
( 2015h also found that stars in the X-shaped bulge do not necessarily stream along 
simple banana or bits. Clearly, more st udies on the orbital structure and vertical reso¬ 
nant heating ('e.g. lOuillen et al.Ll2014l) are desired to make more specific predictions 
for the Milky Way. 
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The detailed kinematics of the near (bright) and far (faint) sides of the X-shape 
may be a useful tool to probe the underlying orbital structure. Several observational 
studies have explore d the X-shape kinematics. With about 300 RCs in (0°,—6°), 
VasQuez et al. ( 2013h found a weak anti-correlation between the longitudinal proper 


motion and radial velocity in both bright and faint RCs. However, they found no 
significant correlation between the latitudinal proper motion and the radial velocity. 
These results were interpreted as possible streamin g motions along the X-shaped 


arms. In bulge fields at b ^ 5°. lPoleski et alJ (l2013h reported the asymmetric mean 


proper motion difference between the near and far sides in both I and b directions; 
this difference is linear for —0.1° < I < 0.5°, but roughly constant for —0.8° < I < 
—0.1°. The linear part was attributed to the streaming motions in the X-shape. 

Numerical simulations can provide comprehensive understanding of the X-shape 
kinematics. The kinemat ic imprints of the X-s hape onto the Galactic bulge observa¬ 
tions were interpreted in [Gardner et al.l (120141) based on their N-body models; they 
found a coherent minimum along 1 = 0° in the mean radial velocity difference be¬ 
tween the near and fa r sides, which is absent in other velocity components and all the 
velocity dispersions. lOin et al.l (l2015h systematically explored the kinematics (both 
the radial velocity and the proper motion) of the X-shape in the SIO model. Along 
I = 0°, the near and far sides of the bar/bulge show excess of approaching and reced¬ 
ing par ticles, reflecting the coherent orbital motion inside the bar structure. Pin et alJ 
( 2015h found only very weak anisotropy in the stellar velocity within the X-shape, 
hinting that the underlying orbita l family of the X-shape may no t be do minated by 
simple banana orbits. Contrary to Poleski et al. ( 2013h . Pin et al. ( 2015h found that 
the proper motion difference between the near/far sides of the X-shape may not be 
used to constrain the bar pattern speed. They also confirmed that the Galactic center 
m ay be located by fittin g the arms of the X-shape, which was originally proposed 
by Gardner et akl ( 2014 ). 

It is important to keep in mind that the true 3D shape of the X-shape should 
not be visualized as simple as a letter “X” with four or eight conspicuous arms 
sticking out. Fig. |6] may give you such an impression because human eyes are eas¬ 
ily biased to pick out the small-scale density enhancement. The true 3D shape 
or structure of iso -density su r faces should really be more like a peanut. This 
is demonstrated in Li & Shen ( 2015h who estimated the 3-D volume density for 
three A-body simulation s with t he adaptiye kernel smoothing technique (ISilyerman . 
1986t Shen & Sellwood . 2004 ): these three bars haye undergone different ampli¬ 
tudes/strengths of buckling instability. Fig. |7] shows clearly that a buckled bar is 
composed of three components with increasing sizes: a central boxy core, a peanut 
bulge and an extended thin bar (Fig. |7]i. The true 3D structure of the X-shape is 
actually more peanut-shaped (Fig. |2]l, but the peanut-shaped bulge can still qualita- 
tiyely reproduce the obseryed bimodal distance distributions that were used to infer 
for the X-shape. Our yisual perception of seeing an “X” is enhanced by the pinched 
concaye shape of the inner peanut structure. 
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Fig. 7 Isodensity surfaces of three W-body bars with different buckling strengths, namely, Model 
1 (top row, strongest buckling), Model 2 (middle row, intermediate buckling, the SIO model) and 
Model 3 (bottom row, negligible buckling). The left and middle columns show the face-on and 
side-on appearance of the bar at large scale (~ 3 kpc). The right column shows the side-on shape 
of the isodensity surfaces at small scale (~ 1 kpc). 


4 More sophisticated chemo-dynamical models of the Galactic 
Bulge 


Despite its simplicity, the self-consistent Milky Way bar/bulge model described in 
Section |2] is successful in explaining many aspects of the Galactic bulge, such as 
the boxy shape, the stellar kinematics, the X-shape, and gives reasonable bar angle, 
length, and the pattern speed. It may also explain, in principle, the vertical metal- 
licity gradient ( Martinez-Valpuesta & Gerhardl2013 and Section |Z2]) . However, in 
order to match the full sophisticated distribution of the stellar populations and chem¬ 
ical composition in the Galactic bulge, more realistic chemo-dynamical simulations 
must eventually supersede the simple A^-body disk simulations (Section lU in the 
future. It is not trivial to model gas properly in simulations, because we need to treat 
many complicated micro-physicical processes of the multi-phase medium with sim- 
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plified prescriptions. In particular, there are still considerable uncertainties in how 
to model the star formation and feedback processes. 

Most chemo-dynamical models are motivated to explain the vertical metallicity 
abundance of the Galactic bulge, which was thought to be the main challenge for 
the secularly evolved bar/bulge model. However, as we discuss in Section 12.21 the 
vertical metallicity gradient may be explained by the simple bar/bulge model as 
well. With additional mass components in the more complicated chemo-dynamical 
simulations, and many more free parameters, the vertical metallicity gradient may 
also be reproduced in these simulations. 


4.1 Two-component disk scenario 


Bekki & Tsuiimoto ( 201 ih reported that their pure thin disk model failed to repro¬ 
duce the vertical met allicity gradient in the Galactic bulge, w hich is opposite to the 
conclusion reached in iMartinez-Valpuesta & Gerhardl(l2013l) . T his difference is due 
mainly to the different initial setup in t he two studies. Unlike Bekki & Tsuiimot^ 
(2011), Martinez-Valpuesta & Gerhard ( 2013 ) did not try to link the initial radial 
metallicity profile to the metallicities of the present-day Galactic disk near the Sun, 
since the buckling instability in the Milky Way must have occurred long ago when 
the outer disk was only incompletely assembled. 

Nev ertheless, in order to explain the observed metallicity gradient Bekki & Tsujimot^ 
( 201 ih proposed a two-component disk scenario where the bulge is formed from 
a disk composed of thin and thick disks. In this scenario, the first thin disk was 
disturbed and heated into a thick disk via a minor merger with a dwarf galaxy 
at early times. A subsequent thin disk is gradually built up with the inner part 
developing a bar structure. Since more metal-poor stars at higher latitude origi¬ 
nate from the already dynamically hotter thick disk, which was not strongly in¬ 
fluenced by vertical mixing of the later bar, they are able to stay in-situ for much 
longer and keep the metallicity low at high vertical distance. Consequently, a ver¬ 
tical metallic ity gradient of the bulg e can be produced. A similar model was also 
suggested in IPi Matteo et al.l (120141) . This scenario may help exp lain the similar¬ 


ity in stellar population s between the bulge and th e thick disk (IMelendez et al 


20081: Alves-Brito et al. . 2010l : Bensbv et al. . 2011). However, it is still not en¬ 
tirely clear whether or not the bulge stars are a distinct po pulation with different 
kinematics and compositions from the thick disk stars (e.g., iLecureur et al.L 12007 


Fulbright, Me William & Rich. 2007: Minniti & Zoccalil 200^ 
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Clump-origin bulges form through mergers of clumps in a primordial galactic disk. 
There have been attempts to link the Milky Way bulge with a clump-origin bulge. 
More discussions of clump-origin bulges may be found in Chapter 6.2 of this book. 

Primordial galaxies are expected to be highly gas-rich in early stage of the disk 
formation, and massive star-forming clumps may form by gravitational instabilities. 
The clump can spiral into the center of the g alaxy rapidly by dynamical friction, th en 
merge to form a central bulge component (Noguchi, 199^ Itnrneli etaL , 2004^ 


u ponent 

^120071: 


.Bournaud, Elmegreen & Elmegree 


Elmegreen, Bournaud & ElmegreenLl2008 


Inoue & Saitohl 120121) . These resu lts are consistent with the clumps in th e chained 
galaxies observed at high redshift ((Elmegreen. Elmegreen & Hirsi 1200411 . 

The reported properties o f the clump-origin bulges still differ consider ably in var¬ 
ious numerical simulations. lElmegreen, Bournaud & Elmegreen! (l2008h found the 
similarity between clump coalescence and major galaxies mergers in terms of or¬ 
bital mixing. In their simulations, these giant star forming clumps migrated towards 
the galaxy center within a few dynamical time scales to form a classical bulge, that 
is thick, rotates slowly and f ollows a Sersic density profile with large index (i.e. 
close to the 7?*/^ law profile). Inoue & Saitohl ( 2012 ) performed high resolution N- 
body/SPH simulations of clump-origin bulges by assuming a collapsing gas sphere 
embedded in an NEW halo. They found that their bulge has many properties indicat¬ 
ing for a pseudobulge, such as a nearly exponential surface density profile, a barred 
boxy shape and strong rotation. They also found that this bulge consists of old and 
metal-rich stars, similar to stars in the Galactic bulge. However, the clump-origin 
bulges made in these studies all failed to produce the definin g characteristic of the 


Milky Way bulge , i.e., the strong cylindrical rotation pattern ((Howard et al.L 12009 


Shen et aUl2010h . The relevance of clump-origin bulges to the Milky Way needs to 
be examined in greater depth in the future. 


4.3 The Milky Way bulge formation in the cosmological setting 

Ideally one would prefer to simulate the full evolution history of Milky Way for¬ 
mation with cosmological initial conditions. Such simulations will simulate the for¬ 
mation of the Galactic bulge in addition to thin and thick disks. The ambitious full 
chemo-dynamical simulations combine 3-D hydrodynamical simulations with the 
calculation of chemical enrichment, and obtain the positions, ages, chemical com¬ 
positions, kinematics of all star particles. Such detailed predictions can be compared 
to observations in large surveys. However, cosmological chemo-dynamical simu¬ 
lations are very computationally expensive, and much realism was still fudged in 
sub-grid physics. 

With improved star formation and feedback models, cosmological disk simula¬ 
tions have made progress in making late-type disk galaxies in recent years. These 
simulations often use strong feedback models to prevent overproduction of stars 
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at early times (see the recent review by Gerhardl 2014 ). Obreia et~ai1 ( 2013h ana¬ 
lyzed and compared the bulges of a sample of L* spiral galaxies in hydrodynamical 
simulations in a cosmological context. They found that the bulges show an early 
starburst-collapse fast phase of mass assembly, followed by a second phase with 
lower star formation rate, driven by disk instabilities and/or minor mergers. They 
suggested that one may associate the old population formed during the first rapid 
phase with a classical bulge, and the young one formed during the slow phase with 
a pseudobulge. The young population is more oblate, generally smaller, more ro- 
tatio nally supported, wit h higher metallicity and less alpha-enhanced than the old 


one. 


Guedes et al. ( 20131) followed the formation and evolution of the pseudobulge 


in the Eris hydrodynamic cosmological simulation with a very high resolution. Their 
pseudobulge was built from the inner disk, and composed of mainly old stars that 
formed in the first step in t he inside-out formation of the baryonic disk. So far these 
cosmological simulations ( Guedes et al. . 20131: Kobavashi . 20141) still do not have 
a boxy/peanut-shaped bulge as in the Milky Way. Cosmological simulations need a 
huge dynamic range to model the bulge and halo simultaneously, since the bulge is 
only ~ 1 kpc in size, and the halo is around ~ 200 kpc. So the spatial resolution re¬ 
quired to resolve structures in the bulge may be too high. There is still a long way to 
go for the current chemo-dynamical simulations to explain the full evolution history 
of the Milky Way’s bulge. 

The ultimate goal of chemo-dynamical simulations is not only to reproduce the 
final results similar to observations, but also to make clear which process is respon¬ 
sible for each particular aspect of observations, which may require a large set of 
parameter studies. Only until we identify the basic factors contributing to the more 
complex features can we claim to understand all the physics involved in the forma¬ 
tion and evolution of the Milky Way. 


5 Other modelling techniques 

5.1 Bar properties constrained from gas kinematics 


The presence of a bar structure in the Mliky W ay was first hinted by the non-circular 
gas kinematics (e.g., de Vaucouleursi 19641) . One may further use the features in 
the asymmetric gas flow pattern to infer the properties of the Galactic bar/bulge. 
Non-circular motions of atomic and molecular gas in the Galaxy are commonly 
presented in the / — v diagram, which shows the distribution of gas emission line 
intensity as a function of Galactic longitude and line-of-si ght velocity since dis¬ 


tances to individual gas clouds are di fficult to measure (e.g.. lBurton & Lis^ll993 


Dame, Hartmann & Thaddeus . 200 ih . The high density features in the I — v diagram 


represent the over-dense regions of gas distribution driven mainly by the large-scale 
non-axisymmetric structures such as the Galactic bar and spiral arms. The features 
in the I — v diagram, due to their unknown distances, must be interpreted through 
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gas dynamical models, and they can provide important constraints on the properties 
of the bar and spiral arms. 

There have been many hydrodynamic models of the gas fl ow in barred potentials 
derived from COBE near-infrared maps or star counts (e.g., Englmuier^GerhM^ 


I 999 I: Bissantz, Englmaier & Gerhard 20031: Rodriguez -Fernan dez & CombesL 2008 ) 


or in self-consistent barred models (e.g., Fux . 1999 : Baba. Saitoh & WadaL 201 (if 
These models have been able to reproduce many of the strong features in the / — v 
diagram, even t hough no m odel was able to provide a good match to all the ob¬ 
served features. Fux ( 1999t) modelled gas dynamics with 3D A^-body - 1 - SPH sim¬ 
ulations. He found that the gas flow driven by a self-consistent bar is asymmetric 
and non-stationary. Snapshots at some specific times can qualitatively reproduce the 
observed / — v diagrams of HI and CO. His model could reproduce the connecting 
arms, which represent the dust lanes as a result of strong bar-driven shocks, 3-kpc 
arm and 135-km s^* arm (both arms are emanating from the ends of the bar). Based 
on his modeling results, he found that a bar angle of 25° ±4°, a bar co-rotation ra¬ 
dius of 4.0-4.5 kpc, and a b ar pa t tern sp eed of ~ 50 km s^'kpc^', consistent with 
hi s stellar modelling resu lts ( Fuxi 1997 ). 

Englmaier & Gerhardl ( 19991) studied gas dynamics in a rotating gravitational po¬ 


tential of the deprojected COBE near-infrared bar and disk. In their models, gas 
formed a pair of shocks at the leading side of the bar and a nuclear ring, typical of the 
gas flow pattern in a barred potential. Many observed gas dynamical features could 
be found in their models, such as the four-armed spiral structure, the 3-kpc arm, the 
terminal velocity curve (i.e., the envelope line of the I — v diagram), tangent points, 
although some features are not exactly identical to those in observations. They pre¬ 
dicted a bar pattern speed of ^ 60 km s~^kpc~* and a bar angle of about 20 — 25°, 
similar to the results from Fux ( 1999 ). Bissantz. Englmaier & Gerhardl ( 20031) sim¬ 
ulated the gas flow with an updated potential further constrained by COBE maps 
and clump giant star counts in several bulge fields. They adopted separate pattern 
speeds for the bar and spiral arms, and suggested that such a configuration would 
help to make spiral arms pass through the bar co-rotation radius where the spiral 
arms dissolve in the single pattern speed simulations. The 3-kpc arm and its far-side 
counterpart in this work cannot be reproduced very well unless they use a massive 
spiral arm potential. They predicted a bar pattern speed of ^ 60 km s^'kpc^' and 
spiral arm pattern speed of ~ 20 km s^'kpc^', similar to previous results. 

The bar pattern speed derived from gas dynamics may have some degeneracy 
with the bar size; a large bar length cou pled with a low er bar pattern speed may 
work quite well to match gas observations. Li et al. ( 2015h modelled gas flow pattern 
for the Milky Way using grid-based hydrodynamical simulations. Their basic bar 
poten tial was from an A^-body model constrained b y the density of bulge red clump 
stars ( Portail et all 2015t Wegg & Gerhard . 2013h . They found that a low pattern 
speed model for the Galactic bar may work, and the best model can give a better fit 
to the / — V diagram than previous high pattern speed simulations. 
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5.2 Orbit-based and particle-based modelling 


Although numerical simulations such as A-body can offer the full evolutionary his¬ 
tory from plausible initial conditions, they also have weaknesses. Numerical sim¬ 
ulations are inflexible in the sense that a lot of trials are required to reproduce the 
desired results by adjusting the initial conditions — this becomes very challenging 
especially for more sophisticated chemo-dynamical simulations which often rely on 
fudged sub-grid physics and are computationally costly, thus limiting the systematic 
exploration of parameter space to match the observational results. 

Other modelling approaches such as the Schwarzschild and made-to-measure 
methods are complementary to A-body simulations. Both methods are good at steer¬ 
ing models to match the desired final results. 

The basic principle of the Schwarzschild orbit-superposition method is given in 
Schwarzschild ( 1979h . One first computes and finds the typical orbital families in 
a potential arising from an assumed triaxial density distribution. The time-averaged 
density along each orbit in a lattice of cells spanning the volume of the model is sim¬ 
ply the total cumulative time spent by that orbit in each cell. Then linear/quadratic 
programming techniques are used to hnd the non-negative weights of each orbit to fit 
the assumed mass distribution and the kinematical data to achieve self-consistency. 
There is no unique solution for the orbital weights, but one may seek to maximize 
some objective function for a balance of minimizations and smoother phase- 


some objective tunction tor a balance ot y minimizations and smoother p hase- 
space distributions (e.g., Richstone & Tremainel 1988t Thomas et al. . 20051) . The 


stability of the constructed model is also not guaranteed, and generally needs to 
be tested with an A-body code. The weighted orbital structure resulting from the 
Schwarzschild modeling can offer important clues to the formation history of sys¬ 
tem. 

Based on the Schwarzschild orbit-superposition method, Zhaol ( 199^ developed 
the first 3D rotating bar model that htted the density prohle of the COBE light dis¬ 
tribution and kinematic data at Baade’s window. His model was constructed with 
485 orbit building blocks, and little stellar kin ematic data were av ailable to ex¬ 
plore the uniqueness of this steady-state model. Hafner et al. ( 2000l) extended the 
classical Schwarzschild technique by combining a distribution function that de¬ 
pends only on classical integrals with orbits that respect non-classical integrals, i.e., 
Schwarzschild’s orbits were used only to represent the difference between the true 
galaxy distribution function and an approximating classical distribution function. 
They used the new method to construct a dynamical model of the inner Galaxy 
with an orbit library that contains about 22,000 regular orbits. For definiteness, they 
assumed a bar angle of 20° and bar pattern speed of 60 km s^'kpc^^ The model 
reproduced the 3D mass density obtained through deprojection of the COBE sur¬ 
face photometry, and the then-available kinematics within the bar corotation radius 
(3. 6 kpc). 

Wang et al. ( 2012 ) extended the Schwarzschild implementation of Zhaol ( 19961) 


and applied it with the extra kinematic constraints from the full BRAVA dataset 
( Kunder et al.L 20121) . Using minimization, their best-fitting Galactic bar model 


has a pattern speed of 60 km s ^kpc *, a disk mass of 10 * ^Mq and a bar angle of 20° 
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out of 36 models varying these parameters. Compared to Zhaol ( 199^ . Wang et ^ 
(120 12h can better reproduce the average radial velocity and the surface brightness 
distribution. However, their model over-predicted the longitudinal proper motions 
compared to the observed values. A^-body tests showed that the model was stable 
only for a short period of 0.5 Gyr. They suspected the ins tability arises beca use no 
self-consistency was imposed for the disk outside 3 kpc. Wang et al. ( 2013h made 
further tests of their implementation using the A^-body bar model in Shen et al] 
(l2010h. wi t h the hope to recover the given bar pattern speed and the bar angle in 


Shen et al. ( 2010l) . They concluded that BRAVA radial velocities alone do not con¬ 


strain well the bar angle and/or the pattern speed using their Schwarzschild im¬ 
plementation, and the observed proper motions may help to reduce the model de¬ 
generacy. Their method appeared to over-fit the BRAVA data points, indicating the 
implemented smoothing of the phase space distributions may need improvement. 

Building a complete and representative orbit library in bars, essential for the 
Schwarzschild method, is non-trivial. It is necessary to explore systematically which 
of the many possible orbit families could contribute to a self-consistent model, and 
to understand how they are affected by properties of the bar potential. For exam¬ 
ple, Schwarzschild modelling may help elucidate the actual orbital compositions 
supporting the prominent X-shaped/peanut-shaped structure in the Ga Hctic bulge . 


There seem to be a large number of irregular/chaotic orbits in bars ('e.g.. lWang et al 


2012 : Patsis & Katsanikas . 2014allb ). We need to have a better idea on the fraction 


of chaotic orbits, classify them into major families, and identify the location of res¬ 
onant orbits. 

Uniqueness of the solution is also an undesired feature in the Schwarzschild mod¬ 
eling of the Galactic bar. There is considerable freedom to reproduce the existing 
data by different sets of orbital weights. More systematic studies are needed to ex¬ 
plore how many different combinations of orbits still reach the same goodness-of-fit 
for a given potential, and how much the properties of the final bar can vary but still 
satisfy the same observational constraints. 


Un like the Schwarzschild method, made-to-measure (M2M) method (ISver & Tremaine , 
1996h slowly adjusts the weights of the particles in an V-body system, instead of the 


orbital weights as in the Schwarzschild method, as particles proceed in their orbits 
until the time-averaged density field and other observables converge to a prescribed 
observational value. The particle weights are adjusted through a weight evolution 
equation according to the mismatch between the model and target observables. In 
the Schwarzschild method orbits are first separately integrated in a fixed potential 
and then superimposed, whereas in the M2M method the two steps are merged — 
orbits are integrated and the weights of particles are adjusted at the same time, thus 
eliminating the need for an orbit library. Also M2M can allow us to dynamically 
adjust the potential while Schwarzschild does not. 

The M2M method was first applied to constr uct a dynamical model for the 
barred bulge and disk model of the Milky Way in IBissantz, Debattista & Gerhard 
( 20041) . Since then the M2M method has been continuously improved in various 
implementations, and has gained growing interests in the dyn amical modelling 


of galaxies, especially spheroidals and early-type galaxies (e.g., Ide Lorenzi et al 
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2007 


2008 200^ Dehnen . 20091: Long & 2010[ 201^ Morganti & Gerhard . 


2012 ; 


Hunt & Kawatal 2013 : Hunt. Kawata & Martel . 1201 A modified j^M2M 


was implemented to improve the algorithm to model both the density and kine¬ 
matic data, account for observational errors and seeing effects, and incorporate a 


maximum-l i kelihood techn i que to account for discrete velocities (Ide Lorenzi et al 


20071 2008). Long & Maol (l2012h a pplied their design and implementation of the 
M 2M method ( Long & M ag 2010l) to 24 SAURON galaxies previously analysed 
by Cappellari et al.l ( 2006 ). and found generally good agreement between M2M and 
Sc hwarzschild metho ds in determining the dynamical mass-to-light ratio. 

(12013h constructed a M2M model of the Galactic bar/bulge con- 


Long et al 


strained by the BRAVA ki nematics ( Kunder et al.L 2012h . They took the A-body 
model in Shen et al. ( 2010l) as initial condition and the target luminosity density. 
They ran a suite of 56 models with different pattern speeds and bar angles in search 
of the best-f itting one. Their b est-fitting model recovered the bar angle and pattern 
speed of the Shen et aD ( 2010l) A-body model, and reproduced both the mean radial 
velocity and radial velocity dispersion of the BR AVA data very wel l. Since they used 
BRAVA results as kinematic constraints and the lShen et al.l ( 201^ model as the tar¬ 
get luminosity density, this work is actually a cross- c heck between the direct A-body 
model ling and the M2M method. Hunt & Kawatal ( 2013 ): Hunt. Kawata & Martell 
( 2013h modelled disk systems with their own design and implementation of the 
M2M method (particle-by-particle M2M method, or PRIMAL). They showed that 
PRIMAL can recover the radial profiles of the surface density, velocity dispersions, 
the rotational velocity of the target disks, the apparent bar structure and the bar 
pattern spee d of the bar. 

Recently, Portail et al. ( 20151) constructed dynamical models of the Galactic 
box/peanut bulge, using the 3D density of red clump giants (IWegg & Gerhardll2013h 
and BRAVA kinemat ics as observational constraint s. They tried to match the data us¬ 
ing their M2M code ( de Lorenzi et all 2007 , 2008 ), starting with A-body models for 
barred disks in different dark matter haloes. In this work, they estimated the total dy¬ 
namical mass (including both stellar and dark matter halo) as 1.84 ± 0.07 x 1O'‘’M0 
inside the rectangular box of ±2.2 x ±1.4 x ±1.2 kpc). They u sed BRAVA kine¬ 


matic a l data to constrain the models, but the proper motion data (IRattenburv et al 


2007 ah are used only as a check of their modelling. Their models tend to be more 


anisotropic than the data. Given the different significance of the disk component 
in their initial conditions, their five models sample a wide range of pattern speeds 
of the final bar structures, ranging from R 1.08 to 1.80, where R ~ RcR/Rbar is 
the dimensionless parameter characterizing the pattern speed of the bar. They got a 
scaled value of 25 to 30 km s^^kpc^*, significantly lower than previous estimates 
with various methods (40 to 60 km s^^kpc^^). If confirmed, this puts the Galactic 
bar among slow rotators (R > 1.5). 

In summary, both the Schwarzschild and M2M methods are designed to steer 
models with pre-determined initial conditions towards prescribed observed results. 
M2M models include aspects of both Schwarzschild methods and regular A-body 
simulations. When the gravitational potential of the target system is held fixed, the 
searching process for a distribution of particle weights is closely related to that for 
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a distribution of orbital weights in Schwarzschilds method to fit the same observa¬ 
tional constraints. Conversely, when the adjustment of particle weights is switched 
off and the p otential is allow ed to evolve, M2M particle codes can reduce to A^-body 
simulations ( Gerhard . 2010l) . On the other hand, both types of modelling only con¬ 
struct self-consistent equilibrium models, but they do not tell us how the Galaxy 
evolved into the current equilibrium configuration from what initial conditions. The 
most important value of these modelling techniques is that they give us predictive 
power and suggest further observational tests to confirm or constrain the structural 
parameters of the bar/bulge. 


6 Summary and future outlook 


Understanding the structure and formation of our Milky Way bulge is nontrivial, 
mostly because of our location in the disk plane and the severe dust extinction in the 
optical band. Near infrared images from the COBE satellite presented the first clear 
evidence of a boxy bulge in the Galaxy. Recent large dedicated surveys have allowed 
a large number of bulge stars to be studied individually in detail. Most bulge stars are 
very old with a wide range of metal abundances, and they formed earlier than most 
disk stars on a rapid formation time scale. Observationa lly, we still need to better un¬ 
derstand the bulge me tallicity componen ts identified by Ness et al. ( 2013ah and their 
kinematic signatures ( Ness et an.l2013bli . Future large suryeys will also need to set¬ 
tle whether or not the recently identified X-s haped structure is populated mainly b y 
metal-rich stars, instead of metal-poor ones (iNataf, Cassisi & AthanassoulaLl2014t) . 

The Galactic bulge contains crucial information about the formation of eyolu- 
tionary history of the Milky Way. To unrayel these clues theoretical modelling of the 
bulge is essential. Here we haye reyiewed recent adyances in modelling the Galactic 
bulge with N-body, chemo-dynamical simulations, and other modelling techniques. 
The main body of the Milky Way bulge appears to be a buckled/thickened bar seen 
somewhat end-on, as hinted from its asymmetric boxy shape. One can construct a 
fully eyolutionary bar model that matches many properties of the Galactic bulge 
reasonably well (Section |2]i. The dynamical eyolution of the bar/bulge was driyen 
mainly by two consecutiye disk instabilities. The bar forms naturally from a cold 
massiye precursor disk yia the well-known bar instability. Shortly after its forma¬ 
tion, the bar suffers from a yigorous buckling instability, and becomes a thickened 
structure that appears boxy or peanut-shaped when seen edge-on. Such a model 
self-consistently eyolyes from plausible simple initial conditions, and is successful 
in explaining many aspects of the Milky Way bulge, such as the excellent match to 
the kinematics of the whole bulge, the X-shaped structure that naturally arises in the 
bar buckling process, reasonable bar angle and other bar parameters consistent with 
independent structural analysis, and the metallicity map. 

This simple model proyides a promising starting point, but there are still many 
open questions to be answered in more sophisticated chemo-dynamical models. For 
example, what is the exact mass fraction of a possible classical bulge hidden un- 
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derneath the dominant boxy bulge; how the fossil record of the early inner Galaxy 
(thick disk, old and younger thin disks) is mapped into the bulge str ucture, i.e., 
how to better understand the bulge metallicity components identified by Ness et aP 
( 2013alf9l in the ARGOS survey, and their correlation with the kinematics; how the 
strongly barred X-shaped structure is populated preferentially by metal-rich stars; is 
there a solid connection betw een the bulge and the thick disk ( Bekki & Tsuiimotol 
201 it Di Matteo et al. . 2014 )? 


Ongoing and upcoming large surveys will undoubtedly shed new light on the 
Milky Way bulge. Gaia will provide accurate parallaxes and prope r motions of 


about 20 million stars along all the lines of sight towards the bulge (IRobin et al 


2005h . Complementary to the Gaia mission, ongoing ground based surveys such 


as APOGEE, VVV, Gaia-ESO, GIBS will also provide us huge amount of high- 
resolution spectroscopic data. With the large influx of data and the improvements 
in theoretical models, we are poised to make greater progress in putting together all 
puzzle pieces of the Milky Way bulge. 
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